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We use a stochastic model and direct numerical simulation to study the impact of turbulence on 
cloud droplet growth by condensation. We show that the variance of the droplet size distribution 
increases in time as f 1 / 2 , with growth rate proportional to the large-to-small turbulent scale separa¬ 
tion and to the turbulence integral scales but independent of the mean turbulent dissipation. Direct 
numerical simulations confirm this result and produce realistically broad droplet size spectra over 
time intervals of 20 minutes, comparable with the time of rain formation. 


Many multi-scale processes- including nutrient forag¬ 
ing of plankton, gas/dust accretion disks in astrophysics 
and spray evaporation and combustion in engines |l|i4]- 
involve the interaction between small particles and trac¬ 
ers transported in a turbulent flow. Here we focus on 
the case of droplet condensation in turbulent warm (i.e. 
ice-free) clouds. Clouds are a leading source of uncer¬ 
tainty in climate modeling Q due to the difficulty of ac¬ 
curately parameterizing the macro-scale effects of micro¬ 
scale physical processes, such as the effect of droplet size 
distribution on precipitation rates and cloud albedos. 

The role of turbulence in cloud microphysics presents 
a range of open questions @-[§|, particularly as a pos¬ 
sible solution for the “bottleneck” problem of rain for¬ 
mation. All cloud droplet populations evolve through 
a sequence of steps: (1) nucleation/activation of cloud 
condensation nuclei, (2) droplet growth by condensation 
and (3) growth to raindrop size by collision and coales¬ 
cence. The passage from (2) to (3) presents a bottle¬ 
neck because collisional growth is only triggered when 
the droplet population acquires a sufficiently broad size 
distribution, but condensational growth is inversely pro¬ 
portional to droplet radius which intrinsically tends to 
narrow the size distribution. Nonetheless, warm clouds 
are routinely observed to precipitate within ~20 minutes 
of formation. Understanding the mechanisms that break 
the condensational bottleneck is a longstanding and still 
unresolved problem in atmospheric physics. 

Turbulence has often been invoked as a key process in 
this context since it can enhance collision rates via iner¬ 
tial clustering 00 and the so called “sling effect” Q 
Turbulence also induces fluctuations in the supersatu¬ 
ration field that can potentially broaden droplet spec¬ 
tra in the condensational stage @]. Early studies using 
analytical models [12), 0] and direct numerical simula¬ 
tions (DNS) 0 generally showed too little broadening as 
compared with observations fl5l |. Later work attempting 
to simulate large-scale turbulence in an 0(100 m) cloud 
113 showed a dramatic broadening of the droplet size 


spectrum but only with ad-hoc assumptions about the 
small-scale supersaturation fluctuations. Lanotte at al. 
19] modelled both small- and large-scale effects on the 


droplet size distribution with simulations of a cloud of 70 
cm and pointed out the importance the Taylor Reynolds 
number, Re\, the non dimensional parameter measur¬ 
ing the large to small scale separation in homogeneous 
isotropic turbulence. In particular, they suggested that 
droplet spectral broadening should scale with Re 

A question that has not been addressed so far is the 
long-term fate of the droplet spectrum: does it reach 
a steady state, or does it continue to evolve? The large 
range of scales involved makes DNS very computationally 
demanding, and all existing simulations consider time 
spans between few seconds and two minutes, well below 
the observed rain formation timescale. 

Here, we approach this question by first deriving an 
analytical expression for the standard deviation of the 
droplet radius squared (droplet surface area) in terms 
of the thermodynamics and turbulence characteristics, 
modeling the fluctuations as stochastic processes. We 
demonstrate that the droplet size distribution increases 
monotonically with time as t 1//2 . We validate this an¬ 
alytical result by extending previous numerical results 
with DNS and large eddy simulations (LES) to timescales 
comparable with those of rain formation, about 20 min¬ 
utes. Our results imply that every warm cloud would 
precipitate if given enough time. The broadening rate is 
determined by the large scale turbulent motions and by 
the positive correlation between droplet surface area and 
local supersaturation. 

Our physical model is analogous to that in 0 , a de¬ 
tailed description can be found in Supplemental Mate¬ 
rials [2cj. We assume homogeneous isotropic turbulence 
obeying the incompressible Navier-Stokes equations 
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where u is the divergence-free fluid velocity, p the pres¬ 
sure, p the air density, f an external forcing to maintain a 
statistically stationary state and v the kinematic viscos¬ 
ity. This approximation is valid for clouds smaller than 
L ss 100 m, which allow spatial inhomogeneity and large 
scale variations of the thermodynamic parameters to be 
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safely neglected. The supersaturation field s is trans¬ 
ported by the fluid according to 


ds s 

— + u • Vs = kV 2 s + A\w - 

at r s 


( 2 ) 


a generalization of the Twomey model [ 27 ] . The dif- 
fusivity of the water vapor in air is denoted by k, 
w is the velocity component in the gravity direction, 
Aiw is a source/sink term of supersaturation resulting 
from the variation in temperature and pressure with 
height. The supersaturation relaxation time r s depends 
on droplet concentration and dimensions | 28 | . t~ x = 
AupwAzAz Ri/y where Ri are the radii of the droplets 
in the volume V, p w the water density, Ai, A 2 and A 3 
thermodynamic coefficients 0- The droplets are as¬ 
sumed to behave as rigid spheres smaller than the Kol¬ 
mogorov turbulent scale, at low mass fraction to neglect 
feedback on the flow. The only forces governing the 
droplet motion are gravity and the Stokes drag (nucle- 
ation/activation is not considered): 
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with x d and v d the droplet position and velocity, u(x d , t) 
the fluid velocity at droplet position, = 2p w R 2 /(9pv) 
the droplet relaxation time, g the gravitational accel¬ 
eration. The supersaturation at the droplet position, 
s(xd,t), determines the droplet evolution via 


^ = 2A 3 s( Xd ,t) (4) 

An exact equation for the average of the droplet radius 
fluctuations is obtained from 0, 
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showing that ((R 2 ) 2 ) increases linearly with time only 
if the correlation ( s'R 2 ) reaches a non-zero statistical 
steady state. 

To quantitatively estimate the droplet growth, we 
adopt a 1-D stochastic model for the velocity fluctu¬ 
ations Wi and the supersaturation field s' of the i-th 
droplet. Such an approach implicitly assumes that the 
small-scale turbulent motions have a negligible influence 
on the macroscopic observables. This assumption is mo¬ 
tivated by previous results revealing the broadening of 
the droplet size distribution with Re\ 0 and fully jus¬ 
tified a posteriori by the numerical simulations presented 
below. 

Homogeneous isotropic turbulence and supersatura¬ 
tion are modelled with two Langevin equations [0 : 


w[(t + dt) = w((t) — dt 
To 
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where v rms is the root mean square of the turbulent ve¬ 
locity fluctuations, £(£) is a zero-mean Gaussian white 
noise, nondimensionalized in order to be ^-correlated in 
time and T 0 the turbulence integral time scale; and 

s f - s 1 '■ 

s'^t + dt) = s'(£) - -f dt -(- A\w[dt - — L T dt+ 
do {T s } 
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for the supersaturation with forcing from the velocity 
field. Here C ws = ( 1 o' s') / ( v rms y/ (s' 2 )) is the normalized 
velocity-supersaturation correlation, (r s ) is the supersat¬ 
uration relaxation time based on the mean droplet radius 
and pit) zero-mean Gaussian white noise, ^-correlated 
in time. Equation J7]) represents a stochastic version of 
the Twomey model. The mean updraft is zero as the 
mean supersaturation (the mean droplet radius does not 
change); entrainment effects [30], collisions and inhomo¬ 
geneities are also not considered to analyze the conser¬ 
vative case when the droplet spectral broadening is only 
induced by supersaturation fluctuations. 

From 0, 0 and 0, assuming (t s ) -C To as in real 
clouds, the fluctuation correlations become 
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Assuming statistically quasi-steady state conditions 
(^ = 0) we find that ( s' 2 ) qs = Afa 2 ms (r s ) 2 and 


(s'R 2 ') qs = 2 A 3 A 2 v 2 rrns (T s ) 2 T 0 = 2A 3 (s , 2 ) qs T 0 (12) 


which give, using 0, 


a R 2 = V%A3AiVrms{Ts){Tot) 1/2 = \j %(s' 2 ) qs A 3 (T 0 t ) 1/2 . 

(13) 


The model shows that the variance of the droplet distri¬ 
bution increases monotonically in a turbulent flow even 
though the supersaturation fluctuations reach a statisti¬ 
cal steady state s qs . The correlation (s'R 2 ) qs , which 
is directly responsible for the radius growth rate, is 
proportional to the level of fluctuations of the super¬ 
saturation field and to the integral scale of the turbu¬ 
lence, see m ■ Expression m can be formulated in 
terms of Kolmo goro v scales since v rm s — Re ^ 2 v q and 
To ~ 0.06Re\r v |29|: 

a R 2 ~ t).7A 3 Aiv l/ 2 (T s )Re\t 1/2 . (14) 
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Note that for t = Tq (short compared with rain forma¬ 
tion time) the lower limit proposed in [19] is recovered, 
ct/j 2 ~ Re 3 / 2 . From (TTlll we note that a R 2 at a fixed 
time depends only on the scale separation represented 
by Re\ and not on the value of the mean dissipation 
inside the clouds. This implies that clouds with differ¬ 
ent dissipation rate and same Reynolds number have an 
identical behavior in terms of droplet growth by conden¬ 
sation. The droplet/turbulence condensation dynamics 
does not depend on the turbulent small scales: the corre¬ 
lation between the supersaturation field and the droplet 
surface area, governing the distribution broadening, is 
determined by the large flow scales. This result is in 
contrast with the belief that the variance of the droplet 
distribution should not grow indefinitely as turbulence 
tends to decorrelate the particle size from the local sat¬ 
uration field [§}. 


To test our predictions, we run simulations by grad¬ 
ually increasing the size of the computational clouds 
from few centimeters to 100 m. The governing equa¬ 
tions (1-4) are solved with a classical pseudo-spectral 
code for the fluid phase coupled with a Lagrangian al¬ 
gorithm for the droplets [3jJ. All cases share the same 
turbulent kinetic energy dissipation e = 10 -3 m 2 s~ 3 , 
a value typically measured in stratocumuli. This cor¬ 
responds to the same small-scale dynamics, with Kol¬ 
mogorov scale 77 = (z/ 3 /e ) 1 / 4 ss 1 mm, Kolmogorov time 
t v = (v/e) 1 / 2 « 0.1 s and velocity v v = tj/t,, « 1 cm/s. 
We examine droplets with 2 different initial radii, 13 /im 
and 5 fj,m, denoted as case 1 and 2, with supersaturation 
relaxation time r s = 2.5 and 7 s, and same concentra¬ 
tion (130 droplets per cm 3 ). The reference temperature 
and pressure are do = 283 K and P = 10 5 Pa, with 
A± = 5 x 10 -4 m _1 , A 2 = 350 m 3 /kg, A 3 = 50 fj,m 2 /s. 
The simulation parameters are reported in tabled] Note 
that simulation DNS El, carrying order 10 9 droplets, is 
to the best of our knowledge the largest Direct Numerical 
Simulation of a warm cloud ever performed. 


Label 

N* L box 

M 

'Urms 

[m/s\ 

T l 

M 

To 

M 

Re\ 

N d 

DNS Al/2 

64 3 008" 

0.035 

2.3 

0.64 

45 

6 x 10 4 

DNS Bl/2 

128 3 0.2 

0.05 

4 

0.95 

95 

9.8 x 10 s 

DNS Cl/2 

256 3 0.4 

0.066 

6 

1.5 

150 

9 x 10 6 

DNS D1 

1024 3 1.5 

0.11 

14 

3 

390 

4.4 x 10 8 

DNS El 

2048 3 3 

0.12 

30 

4 

600 

3. x 10 9 

LES FI 

512 3 100 

0.7 

142 

33 

5000 

1.3 x 10 14 


TABLE I: Parameters of the simulations. The resolution N, 
the cloud size Lbox, the root mean square of the turbulent 
velocity fluctuations v r ms, and Tl = Lbox/v rms approxi- 

mation of the large turbulent scales. To indicates the integral 
time To = (n/2Vr ms ) f[E(k)/k\dk with k the wavenumber 
and E{k) the turbulent kinetic energy spectra j2t|. The total 
number of droplets is indicated by N d . 



FIG. 1: Root mean square of the square droplet radius fluc¬ 
tuations <j R 2 versus time from simulations (symbols) and the 
prediction of the stochastic model (1131) (lines). Inset: correla¬ 
tion {s'R 2 ) from simulations (thin lines) and from Equation 
(1121) (thick lines). 


The time evolution of a R 2 = xj{{R 2 ') 2 ) is shown in 
Figure Q] for all cases investigated. The data confirm the 
predictions from (THU) , i.e. that a R 2 oc f 1 / 2 . 

The correlation ( s'R 2 ) is displayed in the inset of Fig¬ 
ure [T] (thin solid line): In all cases, {s'R 2 ) reaches a 
statistical steady state, fluctuating around the value de¬ 
termining the growth of a R 2. The turbulence creates a 
positive correlation between supersaturation and droplet 
surface area fluctuations that increases as the turbu¬ 
lent scale separation i.e. the cloud size—increases. The 
agreement between the model and the numerical data 
is remarkable for the largest domain sizes where scale 
separation is significant and viscous effects can be ne¬ 
glected. For small Re\, viscous effects are important 
and the stochastic inviscid model overestimates the cor¬ 
rect behavior. For a detailed comparison between DNS 
and stochastic model see Supplemental Materials [Icj. 

To test the model for a larger cloud size, we perform 
a large eddy simulation (LES FI) of a cloud of about 
100 m size. LES can be seen as a good model for our 
problem since it fully resolves the larger flow scale, those 
relevant to droplet condensation/evaporation, as shown 
above. For numerical details see Supplemental Materials 
[20I ]. The Taylor Reynolds number is 5000. The time 
evolution of a R 2 and of {s'R 2 ) are depicted in Fig. [2] 
together with the analytical predictions from (fT21) and 
m The agreement between the two fully validates our 
model. 

To motivate the use of the variance <j 2 r2 to define the 
droplet size distribution we show in Figure [3] that its 
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FIG. 2: Root main square of the square droplet radius fluc¬ 
tuations a R 2 versus time for the LES simulations (symbols) 
and the model (1131) (lines). Inset: correlation (s'R 2 ) from 
the LES (thin line) and the model (1121) (denoted by ”M.” in 
the legend-thick lines). 


probability distribution is Gaussiam in agreement with 
measurements in real stratocumuli [32|. The data in the 
figure refer to the final simulation time (about 20 min- 



FIG. 3: Probability density functions (PDFs) of the square 
radius fluctuations after a simulation time corresponding to 
about 20 minutes (symbols). The lines represent Gaussian 
distributions with the same variance. The inset reports data 
from the LES of a large cloud. 


utes) and are compared to Gaussian distributions of equal 
variance; note that error bars are about the same size as 
the plotting symbols and not visible in the plot. The 
size distribution from the LES of the large cloud (see in¬ 
set) reveals that the Gaussian can be fitted just using the 
value of <j R 2 from the stochastic model also at this higher 
Reynolds number. 

In summary, we have derived an analytical expression 
for the role of turbulence on the dynamics of droplet con¬ 
densation and validated it against large-scale numerical 
simulations. We show that the standard deviation of the 
square droplet radius fluctuations, a R 2 , increases in time 
as f 1 ' 2 ; the growth rate depends linearly on the turbu¬ 
lent scale separation parametrized by Re\. As shown in 
Figure 2, for a cloud with Re\ = 10000—a typical value 
estimated in cumuli with integral scale of 100 m [§}—our 
expression predicts that <r R 2 reaches values in line with 
observations in real clouds (see jUj, their Figure 4) on 
timescales of less than 20 minutes. 

The stochastic approach proposed here may be gener¬ 
alized to consider additional physics and adapted to dif¬ 
ferent micro-scale phenomena in turbulent environments; 
this may also require a numerical solution of the govern¬ 
ing Langevin equations, something still order of magni¬ 
tudes cheaper than a fully-resolved DNS. Indeed our an¬ 
alytical relation predicts numerical results requiring 10 17 
degrees of freedom. From a practical viewpoint, this in¬ 
dicates the promising potential of modeling approaches 
based on PDFs. 

Our results represent a lower limit for the impact of 
turbulence on warm rain formation since real clouds typ¬ 
ically exceed 100 m in scale and are in general non- 
homogeneous, featuring a fluctuating temperature field 
and vigorous entrainment of relatively dry air from out¬ 
side the cloud leading to enhanced supersaturation fluc¬ 
tuations within the cloud. These additional effects would 
lead to even larger values of a R 2 , more than sufficient to 
explain the spectral broadening observed in real clouds. 
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SECTION A. EULERIAN-LAGRANGIAN 
PHYSICAL CLOUD MODEL 

Derivation of the supersaturation equation 

The equation governing the evolution of the supersat¬ 
uration can be derived from first principles, vapor mass 
and energy conservation, and thermodynamics relations. 
Since the supersaturation is a thermodynamic variable 
(as temperature, internal energy and entropy) one can 
write a conservation equation. The original derivation 
by Twomey that can be found in textbooks of atmo- 
sphericphvsics is obtained neglecting the effects of vis¬ 
cosity [l|, 12[ , which we include below in the most general 
case. 

The supersaturation s is defined as s = q v /q vs — 1 
where q v is the vapor mixing ratio and q vs its saturated 
value. The supersaturation rate of change is obtained by 
differentiating its definition: 


the vertical fluid velocity component. The pressure p can 
be obtained from hydrostatic equilibrium: 


Dp _ gp 
Dt R a 6 


(5) 


where R a is the gas costant for dry air. Substituting Eqs. 
© and ((4]) in Eq. © results in: 
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and consequently Eq. © becomes: 

Ds _ J_ D v2 _ Cd_ s + 1 / ee s gL _ ee s g \ 

Dt q vs v Qv q vs q vs \pR v 9 2 c p pR a 0) W 

s + 1 ee s k e L 2g s + 1 ee s L 2 

q vs pR v 0 2 q vs pR v 9 2 c p d 

Assuming s<C 1 and substituting the expression of q vs , 
the previous equation reads: 
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The vapor mixing ratio follows the conservation equation: 


Dq v 
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( 2 ) 


where D v is the thermal diffusivity of water vapor in air 
and Cd is the condensation rate. An evolution equation 
for q vs is obtained from the Clasius-Clapeyron equation 
linking supersaturation pressure e s with temperature 9: 
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where e is the ratio between water and dry air molecular 
weights, p is the thermodynamic pressure, R v is the gas 
costant for water vapor and L is the latent heat of evap¬ 
oration. The temperature 9 obeys energy conservation: 
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where kg is the thermal diffusivity, g the gravity accel¬ 
eration, c p is the specific heat at costant pressure and w 
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The supersaturation relaxation time is 
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where p w is the water density and Ri are the radii of 
the droplets in the volume V. Eq. © can be exactly 
re-written as 
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The only assumption in the model is introduced now, 
assuming the diffusive terms proportional to the laplacian 
of the supersaturation 


fcV 2 s = P^L V 2 gv _ 


ee s 


R v 0 2 


(15) 


obtaining the final equation for the supersaturation field 


—- = A\w —— + fcV 2 s. (16) 

Dt T s 


Given the high Reynolds numbers inside clouds it appears 
reasonable to assume that diffusion and viscous effects 
are much smaller than the convective transport. Indeed, 
we will show in Section B that the diffusive term does 
not affect the supersaturation for large cloud dimensions. 
The term A\w is a source/sink term of supersaturation 
resulting from the variation in temperature and pressure 
with height. 


Model equations and parameters 


The complete Eulerian set of equations for the turbu¬ 
lent velocity field u, the supersaturation s, the temper¬ 
ature fluctuations d, the vapor mixing ratio q v and the 
buoyancy B therefore reads 
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The velocity field u obeys the Navier-Stokes-Boussinesq 
equations, m, where p' is the pressure deviation from 
its hydrostatic value, p is the density of air, f an external 
forcing acting just at the larger turbulent scales, v the 
kinematic viscosity. The external forcing f maintains a 
turbulent field in the steady state regime. 

In the most general cases turbulent clouds are a tran¬ 
sient phenomenon. The transient effects are relevant 
when observing the cloud for its whole life-time, espe¬ 
cially during the cloud generation and dissipation phase. 
The total life-time of a cloud is of the order of hours/days; 
however, if we study droplet condensation for time inter¬ 
vals of the order of minutes, we can consider the sys¬ 
tem at quasi statistical steady-state, as shown by ex¬ 
perimental observations in real clouds as shallow cu¬ 
mulus clouds j3j and marine stratocumulus clouds Q. 


Most of the DNS of clouds available in literature as¬ 
sume steady state turbulence 0-11. The temperature 
follows the advection-diffusion equation (l2Ql) . Equation 
m links the vapor mixing ratio q v , with its saturated 
value q vs = q vs (0,s) function of the temperature and 
the supersaturation field. The momentum equation is 
forced by the buoyancy field B 1 expressed by equation 
(I22p . where Oq is the reference temperature, q v o is the 
vapor mixing ratio at 0o and s = 0, q w is the water mix¬ 
ing ratio and represents the droplet drag on the flow. In 
adiabatic cloud cores, the buoyancy effects are assumed 
to be negligible ficij |. 

The Eulerian fields are coupled with the Lagrangian 
equation for the droplet evolution: 


dvd _ u(x d , t) - v d 
dt r d 


dx d 

dt 


= v d , 


dt 


2 A 3 (9)s{x d ,t) 


ge z, 


(23) 

(24) 

(25) 


with x d and v d the droplet position and velocity, u(x d , t) 
and s(x d , t) the fluid velocity and supersaturation field 
at droplet position, r d = 2p w R 2 /(9pu) the droplet relax¬ 
ation time. The droplets are assumed to be spherical, 
with radius smaller than the Kolmogorov scale in order 
to be considered point-particles. The only forces acting 
on the droplets are the Stokes drag and the gravitational 
force. 


Computational methodology 


The numerical data set described in the manuscript 
has been obtained from Direct Numerical Simulations 
(DNS) using a classical pseudo-spectral method for the 
Eulerian equations coupled with a Lagrangian solver for 
the droplets. For the fluid phase, the Navier-Stokes- 
Boussinesq equations have been integrated in a three- 
periodic domain using a Fourier spectral method with 
the nonlinear terms de-aliased by the 3/2 rule. The solu¬ 
tion is advanced in time using a third-order low-storage 
Runge-Kutta method; specifically, the nonlinear terms 
are computed using an Adam-Bashforth-like approxima¬ 
tion while the diffusive terms are analytically integrated 
11]. The random forcing f is applied isotropically to the 
first shell of wave vectors, with fixed amplitude, constant 
in time and uniformly distributed in phase and directions 
121 ]. For the dispersed phase, we use the point-particle 
approximation as typical droplet sizes are smaller than 
the flow Kolmogorov length. The same Runge-Kutta 
temporal integration used for the carrier phase is adopted 
for the Lagrangian equations. A trilinear interpolation 
scheme is used to compute the flow velocity, supersatu¬ 
ration and temperature at the particle position. The cor¬ 
respondent trilinear extrapolation scheme is adopted to 
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transfer the information from the droplet location (super¬ 
saturation, radius, temperature) into the Eulerian grid. 
The code is fully MPI parallelized to be efficiently run 
on supercomputing infrastructures with a performance 
linearly scaling up to 10 4 cores. 


SECTION B. NUMERICAL MODELS 


Constant coefficient model 


The system of differential and algebraic equations (1- 
14) can be simplified assuming that the coefficients A 1 , 
A 2 , A 3 are only function of the reference temperature 
6 *o, i.e. the temperature fluctuations can be neglected for 
their calculation |7,|9j, and the buoyancy term can be ne¬ 
glected as well 10]. To check and quantify the order of 
magnitude of the temperature fluctuations, we performed 
a preliminary DNS simulation solving the complete sys¬ 
tem (1-14). The simulation matches the parameters de¬ 
scribed in case DNS D 1 of the manuscript. In particular, 
we simulate a cloud of length L\, ox = 1.5 m, velocity 
fluctuations v rms = 0.1 m/s, Taylor Reynolds number 
Re\ = 390 with a resolution of 1024 3 grid points and 
4.4 x 10 8 Lagrangian droplets. From the numerical re¬ 
sults we obtain a negligible value of the temperature root 
mean square 9 rms = 5 x 10~ 3 K. Since, the supersatura¬ 
tion depends on a non-trivial way on the turbulent field 
via the quantities A 2 , A 3 , t s and the Clasius-Clapeyron 
equation, it is relevant to examine the behavior of its 
probability distribution function (pdf). Figure Q] shows 
the pdf of the supersaturations fluctuations for the full set 
of equations (DNS): the behavior is Gaussian and equal 
to the two results from the constant coefficient model 
reported in the manuscript. 


Direct Numerical Simulation (DNS) of the reduced 
model 


As seen in the previous paragraph, we can assume that 
the temperature field is constant in time and space and 
equals to its reference value 9q for the calculation of the 
constants A 3 , A 2 and A 3 while the fluctuations still af¬ 
fect the evaporation-condensation via the supersatura¬ 
tion equation. This has two consequences: the physics 
of the system is simpler to investigate and the numerical 
solver is faster since the equations are more easy to man¬ 
age. Under the simplifying hypothesis the system (1-14) 
reduces to: 



FIG. 1: Probability distribution function of the supersatura¬ 
tion fluctuations extracted by the DNS simulation (symbols) 
and comparison with its Gaussian fit (line). 


V • u 


du 

~dt 

ds 

at 


+ u • Vu 

+ u • Vs 
dv d 


dt 


dx d 


dt 

dt 


0 , 

_W + 

p 


iA7 2 u + f, 


a:V 2 s + A\w — 

u(x d , f) - v d 
Td 


s 

1 

T s 




Vd, 


2 A 3 s(x d ,t) 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 


where A±, A 2 and A 3 are constant that depends on 
the reference temperature 9 3 . The system (15-20) can be 
numerically solved without any further hypothesis, sim¬ 
plifications or models via Direct Numerical Simulations. 
In the manuscript we performed 8 different DNS cases de¬ 
scribed, see Table I, the largest corresponding to a cloud 
of the order of Lb ox = 3 m (DNS El) with a resolution 
of 10 10 degrees of freedom including the Eulerian grid 
points (2048 3 ) and Lagrangian droplets (3 x 10 9 ). 


Large Eddy Simulation (LES) 

Since the largest DNS simulation allows us to describe 
core clouds of the order of 3 meters, we need to employ 
Large Eddy Simulations (LES) to study clouds with size 
of about 100 meters. The LES equations are obtained 
from the DNS system by applying a low-pass filter to the 
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Eulerian fields that remove the smallest scales. In this 
way, just the large scales are evolved while the effects 
of the small scales on the large ones is modeled. The 
generic Eulerian field / can be decomposed in a filtered 
or resolved component / and a residual or subgrid-scale 
component f sgs . Applying the filter operator to the sys¬ 
tem (15-20), we obtain: 


V • u = 0, 

<9u __ Vp 2 — _ 

— — b u • Vu =-b zA7 2 u + f + V • t. 

at p 

ds 

dt 


(32) 

sgs j (33) 


u ■ Vs = kV 2 s -b Aiw —-—b V • J s „ 
r s 

dv d _ u(x d , t ) + u sgs (x d , t) - v d 
dt r d 

dx d 

~dT = Vd ’ 

dR? 

= 2 A 3 (s(x d , t) + s sgs (x d , £)) 


(34) 

- ge z , (35) 

(36) 

(37) 



where f sgs , J sgs , u sgs and s sgs are respectively the sub¬ 
grid stress tensor, the subgrid supersaturation flux, the 
subgrid velocity and the subgrid supersaturation. The 
subgrid terms should be described using suitable closure 
models, which introduce an unavoidable approximation 
of the physics with respect to the DNS. Here, the sub¬ 
grid stress tensor f sgs and the subgrid supersaturation 
flux a re modelled with to the classic Smagorinski 


closure [13} . The subgrid terms appearing in the La- 


grangian droplet equation u sgs and s sgs are set to zero 
so to consider our LES results as a lower limit for the 
droplet size distribution in the cloud. Infact, consider¬ 
ing the temporal evolution of the droplet square radius 
fluctuations: 


dm') 2 ) da 


R 2 


dt 


dt 


4A 3 ((s'R 2 ') + (s' sgs R 2 ')) (38) 



we are implicitly assuming that the subgrid correlations 
( s' sgs R 2 ) are much smaller than the resolved scale corre¬ 
lations ( s'R 2 ). This is justified by the DNS results rep¬ 
resented in the inset of figure 1 of the main manuscript 
where the correlation (s'R 2 ) is displayed. The figure 
shows that the large DNS gives values of this correlation 
two order of magnitude larger than that from the smaller 
simulations. 

Another approximation introduced in the LES is re¬ 
lated to the number of droplets tracked in the simula¬ 
tion. Theoretically, one should evolve order 10 14 droplets 
which is unfeasible. We therefore use the method of 
droplet renormalization described in Q. We validated 
the LES by comparing the results with our finest DNS 
simulation, DNS El. Different LES simulations have 
been performed with different resolutions down to 32 3 
grid points. The root mean square of the square droplet 
radius fluctuations (top panel) and of the supersaturation 


FIG. 2: Top panel: comparison between the root mean square 
of the square droplet radius fluctuations <j R 2 versus time for 
the LES (symbols) and the DNS simulations. Bottom panel: 
comparison between the root mean square of the supersatu¬ 
ration fluctuations versus time for the LES (symbols) and the 
DNS simulations. 


fluctuations (bottom panel) are shown in figure [2] for the 
DNS and LES. The LES slightly underestimates the ob¬ 
servable a 2 R 2 and s rms as expected from equation (1381) . 
The agreement between LES and DNS is excellent consid¬ 
ering that we used 10 5 degrees of freedom against 10 10 for 
the DNS. This implies that just solving the larger scales 
is enough to provide a good prediction of the droplet size 
spectra variance. In conclusion, the LES simulation is 
able to give a lower limit of the droplet spectra broaden¬ 
ing at a moderate computational cost. 
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FIG. 3: Top panel: Relative error pertaining the quasi steady- 
state supersaturation fluctuation variance in the stochastic 
model (SM) and DNS. Bottom panel: ratio between the 
droplet sink term and the source term in the supersaturation 
variance equation (28). 


Stochastic model (SM) 

The other approach we used in our work is the stochas¬ 
tic model described by the equations (6-7) in the main 
manuscript. We show in the main text that the stochastic 
model tends to overestimate the droplet-related quanti¬ 
ties and that the differences become smaller by increasing 
the ratio between the large and the small scales in the nu¬ 
merical simulations. A systematic study of the difference 
between the stochastic model and the DNS is presented 
in this section. 


The relative error pertaining the quasi steady-state 
value of the supersaturation fluctuation variance is re¬ 
ported in the top panel of figure [3] versus the cloud size 
domain (the size of the computational domain and the 
large-to-small scale ratio). A large error is evident for 
clouds of smaller size and the differences are monotoni- 
cally decreasing when increasing the large-to-small scale 
ratio (i.e. the Reynolds number). For clouds that are of 
the order of 100 meters, our prediction has a margin of 
error lower than 20 %, an acceptable value for statistical 
observables that DNS cannot calculate with the current 
computational resources. 

The error in the stochastic model can be explained by 
viscous effects in the supersaturation fluctuations. Mul¬ 
tiplying equation (17) with the supersaturation fluctua¬ 
tions and averaging, we obtain an equation for the su¬ 
persaturation variance: 

+ V • <uV 2 /2) = -k ((Vs') 2 ) + A^w's') - <— ) 

Ot T s 

(39) 

where the left hand side is zero at statistically steady 
state, and a balance between three terms holds: viscous 
dissipation ac((Vs / ) 2 ), production Ai(w's') and droplet 
sink ( s' 2 /t s ). In the case of the classic inviscid Twomey 
equation (employed also in the stochastic model) the bal¬ 
ance is only between production and the droplet sink 
terms, and viscous effects are neglected. To quantify the 
differences between the inviscid and viscid Twomey equa¬ 
tion, we display the ratio between the droplet sink and 
the production term as extracted from the DNS in the 
bottom panel of figure[3l where a value of this ratio equal 
to 1 corresponds to an inviscid behavior as assumed in 
the Towmey model. For smaller clouds, the ratio is about 
0.5 so that just half of the supersaturation fluctuations 
induce the droplet growth while the remaining is diffused. 
This explain why the stochastic model, neglecting diffu¬ 
sion, is not able to predict the DNS results. By increasing 
the Reynolds number and the size of the cloud, the ratio 
(s ' 2 / t s ) / A\{w's') tends to 1, meaning that the viscous 
effects become more and more negligible (they vanish in 
the limit of infinite Reynolds number). Therefore, the 
stochastic model effectively predicts an upper limit for 
the droplet spectral broadening at a zero computational 
cost. In conclusion, the DNS results will always lie be¬ 
tween the LES estimation (lower limit) and the stochastic 
model prediction (upper limit). 
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